Determining the density of states for classical statistical models: A random walk 

algorithm to produce a flat histogram 
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We describe an efficient Monte Carlo algorithm using a random walk in energy space to obtain a 
very accurate estimate of the density of states for classical statistical models. The density of states 
is modified at each step when the energy level is visited to produce a flat histogram. By carefully 
controlling the modification factor, we allow the density of states to converge to the true value very 
quickly, even for large systems. From the density of states at the end of the random walk, we can 
estimate thermodynamic quantities such as internal energy and specific heat capacity by calculat- 
ing canonical averages at essentially any temperature. Using this method, we not only can avoid 
repeating simulations at multiple temperatures, but can also estimate the Gibbs free energy and 
entropy, quantities which are not directly accessible by conventional Monte Carlo simulations. This 
algorithm is especially useful for complex systems with a rough landscape since all possible energy 
levels are visited with the same probability. As with the multicanonical Monte Carlo technique, our 
method overcomes the tunneling barrier between coexisting phases at first-order phase transitions. 
In this paper, we apply our algorithm to both 1st and 2nd order phase transitions to demonstrate 
its efficiency and accuracy. We obtained direct simulational estimates for the density of states for 
two-dimensional ten-state Potts models on lattices up to 200 x 200 and Ising models on lattices up 
to 256 X 256. Our simulational results are compared to both exact solutions and existing numerical 
data obtained using other methods. Applying this approach to a 3D ±J spin glass model we esti- 
mate the internal energy and entropy at zero temperature; and, using a two-dimensional random 
walk in energy and order-parameter space, we obtain the (rough) canonical distribution and energy 
landscape in order-parameter space. Preliminary data suggest that the glass transition tempera- 
ture is about 1.2 and that better estimates can be obtained with more extensive application of the 
method. This simulational method is not restricted to energy space and can be used to calculate 
the density of states for any parameter by a random walk in the corresponding space. 
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I. INTRODUCTION 

Computer simulation now plays a major role in sta- 
tistical physics [|l|, particularly for the study of phase 
transitions and critical phenomena. The standard Monte 
Carlo (MC) method has been the Metropolis importance 
sampling algorithm |^, but more recently new, efficient 
algorithms have begun to play a role in allowing simu- 
lation to achieve the resolution which is needed to accu- 
rately locate and characterize phase transitions |Q]. For 
example, cluster flip algorithms, beginning with the sem- 
inal work of Swendsen and Wang ||3|] and extended by 
Wolff [Q , have been used to reduce critical slowing down 
near 2nd order transitions. Similarly, the multicanoni- 
cal ensemble method 1q] was introduced to overcome the 
tunneling barrier between coexisting phases at 1st or- 
der transitions and has general utility for systems with 
a rough energy landscape |^-|§]. In both situations, his- 
togram re-weighting techniques |^ can be applied in the 
analysis to increase the amount of information that can 
be gleaned from simulational data, but the applicability 
of re- weighting is severely limited in large systems by the 
statistical quality of the "wings" of the histogram. This 
latter effect is quite important in systems with compet- 
ing interactions for which short range order effects might 



occur over very broad temperature ranges or even give 
rise to frustration that produces a very complicated en- 
ergy landscape and limits the efficiency of "standard" 
methods. 

One of the most important quantities in statistical 
physics is the density of states g{E), i.e. the number of 
all possible states (or configurations) for an energy level 
E of the system, but direct estimation of this quantity 
has not been the goal of simulations. Instead, most con- 
ventional Monte Carlo algorithms |l| such as Metropolis 
importance sampling, Swendsen- Wang cluster flipping, 
etc. generate a canonical distribution g{E)e^^^^'^'^ at 
a given temperature. Such distributions are so narrow 
that, with conventional Monte Carlo simulations, multi- 
ple runs are required if we want to know thermodynamic 
quantities over a significant range of temperatures. In 
the canonical distribution, the density of states does not 
depend on the temperature at all. If we can estimate the 
density of states g{E) with high accuracy for all energies, 
we can then construct canonical distributions at essen- 
tially any temperature. For a given model in statistical 
physics, once the density of states is known we can calcu- 
late the partition function as Z = '^g{E)e~^^ , and the 

E 

model is essentially "solved" since most thermodynamic 
quantities can be calculated from it. Though computer 



simulation is already a very powerful method in statisti- 
cal physics [Q , it seems that there is no efficient algorithm 
to calculate the density of states very accurately for large 
systems. Even for exactly solvable models such as the 2- 
dim Ising model, g(E) is impossible to calculate exactly 
for a large system [Q . 

The multicanonical ensemble method j^jl^ proposed 
by Berg et at has been proved to be very efficient 
in studying first-order phase transitions where simple 
canonical simulations have difficulty overcoming the tun- 
neling barrier between coexisting phases at the transition 
temperature [p|jqjl^-[l^ . The method also has been suc- 
cessfully applied to some complex systems, such as spin 
glass models |0,|lJ-|l8[ and protein folding problems [||, 
for which the energy landscape is very rough and the 
conventional canonical Monte Carlo simulation gets eas- 
ily trapped in local nfinima at low temperature. In the 
multicanonical method, we have to estimate the den- 
sity of states g{E) first, then perform a random walk 
with a probability of -tW to make the histogram ffat in 
the desired region in the phase space, such as between 
two peaks of the canonical distribution at the first-order 
transition temperature. In a multicanonical simulation, 
the density of states need not necessarily be very ac- 
curate, as long as the simulation generates a relatively 
flat histogram and overcomes the tunneling barrier in en- 
ergy space. This is because the subsequent re- weighting 
P,pl| , does not depend on the accuracy of the density of 
the states as long as the histogram can cover all impor- 
tant energy levels with sufficient statistics. (If the den- 
sity of states could be calculated very accurately, then 
the problem would have been solved in the first place 
and we needn't perform any further simulation such as 
with the multicanonical simulational method.) Berg et 
al. proposed a recursive method to calculate the density 
of states by accumulating histogram entries and estimat- 
ing the density of states iteratively with the histogram 
data [0,^- Their method works well for small systems 
since the number of all possible states is small; however, 
for large systems the number of possible states increases 
exponentially with the size of a system. A simple iter- 
ation method to construct a histogram gets trapped in 
a narrow energy range, and it needs an extremely long 
time and large number of iterations to get an accurate 
estimate for the density of states for the entire energy 
space, even for a small system. It is thus not practical 
to calculate the density of states for large systems with 
this approach. For a simple model such as the Potts 
model. Berg and Neuhaus used finite-size scaling theory 
to "guess" the density of states up to 100 x 100 from 
the simulational results of the small systems |j]. Such 
results are not as reliable as direct simulations, moreover 
for complex systems such as spin glass models, we can 
not simply apply such finite-size scaling at all. Berg and 
Celik applied their multicanonical method to a 2D spin- 
glass model, but they could test their method only up to 
a 48 X 48 lattice 0. Recently Berg and Janke proposed 



a similar method (multioverlap simulational) for the 3D 
Ising spin-glass model, they obtained some reliable re- 
sults for systems as large as 12 x 12 x 12 p5[ . 

Lee [ [20[ independently proposed the entropic sampling 
method, which is basically equivalent to multicanonical 
ensemble sampling. He used an iteration process to cal- 
culate the microcanonical entropy at E which is defined 
by S{E) = ln[g(_E)] where g{E) is the density of states. 
He also applied his method to the 2D ten-state {Q — 10) 
Potts model and the 3D Ising model; however, just as 
for other simple iteration methods, it works well only 
for small systems. He obtained a good result with his 
method for the 24 x 24 2-dim Q = 10 Potts model and 
the 4x4x4 3-dim I sing model. 

Oliveira et al. ||2l|-p3|l proposed the broad histogram 
method with which they calculated the density of states 
by estimating the probabilities of possible transitions be- 
tween all possible states of a random walk in energy 
space. Using simple canonical average formulae in statis- 
tical physics, they then calculated thermodynamic quan- 
tities for any temperature. Though the authors believed 
that the broad histogram relation is exact, their simu- 
lational results have systematic errors even for the Ising 
model on a lattice as small as 32 x 32 in references pl|,p5| ■ 
They believed that the error was due to the particular dy- 
namics adopted within the broad histogram method P5| , 
but other work P9,p^ argued that it violates the de- 
tailed balance condition. The algorithm was corrected in 
reference [g8[ and an approach to the broad-histogram 
method was proposed to calculate the density of states 
based on the number of potential moves during the ran- 
dom walk in energy space [24). They obtained more ac- 
curate results for a 32 x 32 Ising model than the broad 
histogram did; but this method also suffers from the sys- 
tematic errors and substantial deviations when system 
becomes larger than 32 x 32 ||H. 

It is thus an extremely difficult task to calculate den- 
sity of states directly with high accuracy for large sys- 
tems. All methods based on accumulation of histogram 
entries, such as the histogram method of Ferrenberg and 
Swendsen [^1, Lee's version of multicanonical method (en- 
tropic sampling) J20| , broad histogram method pl| , p^ , p9[ 
and ffat histogram method ||2J] have the problem of scal- 
ability for large systems. These methods suffer from sys- 
tematic errors when systems are large, so we still need a 
superior algorithm to calculate the density of states for 
large systems. 

Very recently, we introduced a new, general, efficient 
Monte Carlo algorithm that offers substantial advantages 
over existing approaches ||3y]. In this paper, we will ex- 
plain the algorithm in detail, including our implementa- 
tion, and describe its application not only to 1st and 2nd 
order phase transitions, but also to a 3D spin glass model 
that has a rough energy landscape. 

Unlike conventional Monte Carlo methods that directly 
generate a canonical distribution at a given temperature 
g{E)e~^/^^'^ , our approach is to directly estimate the 
density of states g{E) accurately via a random walk that 



produces a flat histogram in energy space. We modify 
density of states at each step of the random walk, and by 
carefully controlling the modification factor we can ob- 
tain a density of states that converges to the real value 
very quickly even for large systems. The resultant den- 
sity of states is accurate enough to calculate thermody- 
namic quantities by applying canonical average formulas 
in statistical physics. 

The remainder of this paper is arranged as follows. In 
section II, we present our general algorithm in detail. In 
section III, we apply our method to the 2D Q = 10 Potts 
model which has a first-order phase transition. In section 
IV, we apply our method to a model with a second-order 
phase transition to test the accuracy of the algorithm. 
In section V, we consider the 3D ± J spin glass model, a 
system with rough landscapes. Discussion and the con- 
clusion are presented in section VI. 



II. A GENERAL AND EFFICIENT ALGORITHM 

TO ESTIMATE THE DENSITY OF STATES 

WITH A FLAT HISTOGRAM 

Our algorithm is based on the observation that if we 
perform a random walk in energy space by fiipping spins 
randomly for a spin system, and the probability to visit 
a given energy level E is proportional to the reciprocal of 
the density of states -rgv, then a fiat histogram is gen- 
erated for the energy distribution. This is accomplished 
by modifying the estimated density of states in a system- 
atic way to produce a "flat" histogram over the allowed 
range of energy and simultaneously making the density 
of states converge to the true value. We modify the den- 
sity of states constantly during each step of the random 
walk and use the updated density of states to perform a 
further random walk in energy space. The modification 
factor of the density of states is controlled carefully, and 
at the end of simulation the modification factor should 
be very close to 1 which is the ideal case of the random 
walk with the true density of states. 

At the very beginning of our simulation, the density of 
states is a priori unknown, so we simply set all entries to 
g{E) — 1 for all possible energies E. Then we begin our 
random walk in energy space by fiipping spins randomly 
and the probability at a given energy level is proportional 
to -tW. In general, if Ei and E2 are energies before and 
after a spin is fiipped, the transition probability from 
energy level Ei to E2 is: 



p{Ei -^ E2) = min 



9{E2y 



(1) 



Each time an energy level E is visited, we modify the 
existing density of states by a modification factor / > 1 , 
i.e. g{E) -^ g{E) * f. (In practice, we use the formula 
ln[g(i<^)] -^ ln[5(i?)] -|- ln(/) in order to fit all possible 
g{E) into double precision numbers for the systems we 
will discuss in this paper.) If the random walk rejects 



a possible move and stays at the same energy level, we 
also modify the existing density of states with the same 
modification factor. Throughout this study we have used 
an initial modification factor of / = /o = e^ ~ 2.71828... 
which allows us to reach all possible energy levels very 
quickly even for a very large system. If /o is too small, 
the random walk will spend an extremely long time to 
reach all possible energies. However, too large a choice 
of /o will lead to large statistical errors. In our simula- 
tions, the histograms are generally checked about each 
10000 MC sweeps. A reasonable choice is to make /o°"°" 
have the same order of magnitude as the total number of 
states {Q^ for a Potts model). During the random walk, 
we also accumulate the histogram H{E) (the number of 
visits at each energy level E) in the energy space. When 
the histogram is "fiat" in the energy range of the ran- 
dom walk, we know that the density of states converges 
to the true value with an accuracy proportional to that 
modification factor ln(/). Then we reduce the modifica- 
tion factor to a finer one using a function like /i = ^/Jo, 
reset the histogram, and begin the next level random 
walk during which we modify the density of states with 
a finer modification factor /i during each step. We con- 
tinue doing so until the histogram is "fiat" again and 
then reduce the modification factor /i+i = ^/Jl ^^'^ 
restart. We stop the random walk when the modifica- 
tion factor is smaller than a predefined value (such as 
/final = exp(10"^) ~ 1.00000001). It is very clear that 
the modification factor acts as a most important control 
parameter for the accuracy of the density of states during 
the simulation and also determines how many MC sweeps 
are necessary for the whole simulation. The accuracy of 
the density of states depends on not only /final but also 
many other factors, such as the complexity and size of the 
system, criterion of the flat histogram and other details 
of the implementation of the algorithm. 

It is impossible to obtain a perfectly flat histogram 
and the phrase "flat histogram" in this paper means that 
histogram H{E) for all possible E is not less than x% 
of the average histogram {H{E)), where x% is chosen 
according to the size and complexity of system and the 
desired accuracy of the density of states. For the L = 32, 
2D Ising model with only nearest-neighbor couplings, this 
percentage can be chosen as high as 95%, but for large 
systems the criterion for "flatness" may never be satisfied 
if we choose too high a percentage and the program may 
run forever. 

One essential constraint on the implementation of the 
algorithm is that the density of states during the random 
walk converges to the true value. The algorithm proposed 
in this paper has this property. The accuracy of the 
density of states is proportional to ln(/) at that iteration; 
however, ln(/fi„ai) can not be chosen arbitrary small or 
the modified ln[(7(i?)] will not differ from the unmodified 
one to within the number of digits in the double precision 
numbers used in the calculation. If this happens, the 
algorithm no longer converges to the true value, and the 
program may run forever. Even if /snai is within range 



but too small, the calculation might take excessively long 
to finish. 

We have chosen to reduce the modification factor by a 
square root function, and / approaches 1 as the number 
of iterations approaches infinity. In fact, any function 
may be used as long as it decreases / monotonically to 1. 

1 In 

A simple and efhcient formula is /i+i — j^ , where n > 
1 . The value of n can be chosen according to the available 
CPU time and expected accuracy of the simulation. For 
the systems that we have studied the choice of n = 2 
yielded good accuracy in a relatively short time, even for 
large systems. When the modification factor is almost 
1 and the random walk generates a uniform distribution 
in energy space, the density of states should converge to 
the true value for the system. 

The method can be further enhanced by performing 
multiple random walks, each for a different range of en- 
ergy, either serially or in parallel fashion. We restrict 
the random walk to remain in the range by rejecting any 
move out of that range. The resultant pieces of the den- 
sity of states can then be joined together and used to 
produce canonical averages for the calculation of ther- 
modynamic quantities at essentially any temperature. 

We should point out here that during the random walk 
(especially at the early stage of iteration) in the energy 
space, the algorithm does not satisfy the detailed bal- 
ance condition exactly, since the density of states is mod- 
ified constantly during the random walk. However, after 
many iterations, the density of states converges to the 
true value very quickly as the modification factor ap- 
proaches 1. li p{Ei — > E2) is the transition probability 
from the energy level Ei to level E2, from equation (|^), 
the ratio of the transition probabilities from Ei to E2 
and from E2 to Ei can be calculated very easily as: 



pjEi ^E2) _ gjEi) 
p{E2 -> El) g{E2) 



(2) 



where g{E) is the density of states. In another words, 
our random walk algorithm satisfies the detailed balance 
condition: 



9{Ei) 



p{Ei ^ E2) = 



1 



giE2 



ME2 ^ El) 



(3) 



where ,\, , is the probability at the energy level Ex and 
p{Ei -^ E2) is the transition probability from Ei to E2 
for the random walk. We conclude that the detailed bal- 
ance condition is satisfied with accuracy proportional to 
the modification factor ln(/). 

Almost all recursive methods update the den- 
sity of states by using the histogram data directly 
only after enough histogram entries are accumu- 
lated ||,|,^3[ll5|-|lf,||,|l|-|3[. Because of the exponential 
growth of the density of states in energy space, this pro- 
cess is not efficient because the histogram is accumulated 
linearly. In our algorithm, we modify the density of states 
at each step of the random walk, and this allows us to 



approach the true density of states much faster than con- 
ventional methods especially for large systems. (We also 
accumulate histogram entries during the random walk, 
but we only use it to check whether the histogram is flat 
enough to go to the next level random walk with a finer 
modification factor.) 

We should point out here that the total number of 
configurations increases exponentially with the size of 
the system; however, the total number of possible en- 
ergy levels increases linearly with the size of system. It 
is thus easy to calculate the density of states with a ran- 
dom walk in energy space for a large system. In this 
paper for an example, we consider the Potts model on 
a, L X L lattice with nearest- neighbor interactions [ p4[ . 
For Q > 3, the number of possible energy levels is about 
2N, where N = L^ is the total number of the lattice site. 
However, the average number of possible states (or con- 

qN 

figurations) on each energy level is as large as ^^ , where 
Q is the number of possible states of a Potts spin and 
Q'^ is the total number of possible configurations of the 
system. This is the reason why most models in statistical 
physics are well defined, but we can not simply use our 
computers to realize all possible states to calculate any 
thermodynamic quantities, this is also the reason why 
efficient and fast simulational algorithms are required in 
the numerical investigations. 

By the end of simulation, we only obtain relative den- 
sity, since the density of states can be modified at each 
time it is visited. We can apply the condition that the to- 
tal number of possible states for the Q state Potts model 
i^ Tli9{-E) = Q^ or the number of ground state is Q to 

E 

get the absolute density of states. 



III. APPLICATION TO A FIRST ORDER PHASE 
TRANSITION 

A. Potts model and its canonical distribution 

In this section, we apply our algorithm to a model 
with a first-order phase transition |g^,Q. We choose 
the 2D ten state Potts model since it serves as an ideal 
laboratory for temperature-driven first-order phase tran- 
sitions. Since some exact solutions and extensive simula- 
tional data are available, we have ample opportunity to 
compare our results with other values. 

We consider the 2-dimensional (5 = 10 Potts model on 
L X L square lattice with nearest-neighbor interactions 
and periodic boundary conditions. The Hamiltonian for 
this model can be written as: 



W = - ^ Siq^,qJ] 



(4) 



<y> 



and q = 1,2, ...Q. During the simulation, we select lattice 
sites randomly and choose integers between [1 : Q] ran- 
domly for new Potts spin values. The modification factor 



fi changes from fo — e^ — 2.71828 at the very beginning 
to /final = exp(10"^) ~ 1.00000001 by the end of the ran- 
dom walk. At the end of the simulations, our algorithm 
only provides a relative density of states for different en- 
ergies, so to extract the correct density of states, we can 
either use the fact that the total number of possible states 
is Q'^ or that the number of ground states is Q, where 
N — L^ is the total number of lattice sites. (Actually we 
can use one of these two conditions to get the absolute 
density of states, and use the other condition to check 
the accuracy of our result.) To guarantee the accuracy 
of thermodynamic quantities at low temperatures in fur- 
ther calculations, in this paper we use the condition that 
the number of the ground states is Q to normalize the 
density of states. The densities of states for 100 x 100, 
150 X 150 and 200 x 200 are shown in Fig. |l|. It is very 
clear from the figure that the maximum density of states 
for L — 200 is very close to 10"'°"°'' which is actually 
about 5.75 x lO'^^^^'^ from our simulational data. 

Conventional Monte Carlo simulation (such as 
Metropolis sampling ||l|,|[) realizes a canonical distribu- 
tion P{E^ T) by generating a random walk Markov chain 
at a given temperature: 



P{E,T)^g{E)e'^'^^'^ 



(5) 



From the simulational result for the density of states 
g{E), we can calculate the canonical distribution by the 
above formula at essentially any temperature without 
performing multiple simulations. In Fig. y(a), we show 
the resultant double peaked canonical distribution [ p6[ , 
at the transition temperature T^ for the first-order tran- 
sition of the Q = 10 Potts model. The "transition tem- 
peratures" Tc(L) are 0.70171 for L = 60, 0.70143 for 
L = 80 and 0.70135 for L — 100 which are determined 
by the temperatures where the double peaks are of the 
same height. Note that the peaks of the distributions are 
normalized to 1 in this figure. The valley between two 
peaks is quite deep. e.g. is 7 x 10~^ for L = 100. The 
latent heat for this temperature driven first-order phase 
transition can be estimated from the energy difference 
between the double peaks. Our results for the locations 
of the peaks are listed in the Table 1. They are consistent 
with the results obtained by multicanonical method |0] 
and multibondic cluster algorithm |^ for those lattice 
sizes for which these other methods are able to gener- 
ate estimates. As the table shows, our method produces 
results for substantially larger systems than have been 
studied by these other approaches. 

Because of the double peak structure at a first-order 
phase transition, conventional Monte Carlo simulations 
are not efficient since an extremely long time is required 
for the system to travel from one peak to the other in en- 
ergy space. With the algorithm proposed in this paper, 
all possible energy levels are visited with equal proba- 
bility, so it overcomes the tunneling barrier between the 
coexisting phases in the conventional Monte Carlo sim- 
ulations. The histograms for L — 60, 80 and 100 are 



shown in the inset of the Fig. ||(a) and are very flat. 
The histogram in the figure is the overall histogram de- 
fined by the total number of visits to each energy level 
for the random walk. Here, too, we choose the ini- 
tial modification factor /o = e^, and the final one as 
exp(10-**) ~ 1.00000001; and the total number of iter- 
ations is 27. In our simulation, we do not set a prede- 
termined number of MC sweeps for each iteration but 
rather give the criterion which the program checks pe- 
riodically. Generally, the number of MC sweeps needed 
to satisfy the criterion increases as we reduce the modi- 
fication factor to a finer one, but we cannot predict the 
exact number MC sweeps needed for each iteration be- 
fore the simulation. We believe that it is preferable to 
allow the program to decide how great a simulational ef- 
fort is needed for a given modification factor fi . This also 
guarantees a sufficiently flat histogram resulting from a 
random walk which in turn determines the accuracy of 
the density of states at the end of the simulation. We 
nonetheless need to perform some test runs to make sure 
that the program will flnish within a given time. The en- 
tire simulational effort used was about 1.6 x 10^ visits per 
energy level (or 3.2 x 10^ MC sweeps) for L = 60, 2.2 x lO"^ 
visits for L = 80 and 3.3 x 10^ visits for L = 100. With 
the program we implemented, the simulation for L = 100 
can be completed within two weeks in a single 600 MHz 
Pentium HI processor. 

To speed up the simulation, we needn't constrain our- 
selves to performing a single random walk over the entire 
energy range with high accuracy. If we are only inter- 
ested in a speciflc temperature range, such as near Tc , 
we could first perform a low precision unrestricted ran- 
dom walk, i.e. over all energies, to estimate the required 
energy range, and then carry out a very accurate ran- 
dom walk for the corresponding energy region. The inset 
of Fig. ||(a) only shows the histograms for the extensive 
random walks in the energy range between E/N = —1.90 
and —0.6. If we need to know the density of states more 
accurately for some energies, we also can perform sepa- 
rate simulations, one for low energy levels, one for high 
energy levels, the other for middle energy which includes 
double peaks of the canonical distribution at Tc. This 
scheme not only speeds up the simulation, but also in- 
creases the probability of accessing the energy levels for 
which both maximum and minimum values of the distri- 
butions occur by performing the random walk in a rela- 
tively small energy range. If we perform single random 
walk over all possible energies, it will take a long time to 
generate rare spin configurations. Such rare energy levels 
include the ground energy level or low energy levels with 
only few spins with different values and high energy lev- 
els where all, or most, adjacent Potts spins have different 
values. 

With the algorithm in this paper, if the system is not 
larger than 100 x 100, the random walk on important en- 
ergy regions (such as that which includes the two peaks 
of the canonical distribution at Tc) can be carried out 
with a single processor and will give an accurate den- 



sity of states within about lO'' visits per energy level. 
The results shown in Fig. |^(a) were obtained using a sin- 
gle processor. However, for a larger system, we can use 
a parallelized algorithm by performing random walks in 
different energy regions, each using a different processor. 
We have implemented this approach using PVM with a 
simple master-slave model and can then obtain an ac- 
curate estimate for the density of states with relatively 
short runs on each processor. The densities of states for 
150 X 150 and 200 x 200, shown in Fig. |l|, were obtained by 
joining together the estimates obtained from 21 indepen- 
dent random walks, each constrained within a different 
regions of energy. The histograms from the individual 
random walks are shown in the inset of Fig. ||(b) both 
for 150 X 150 and 200 x 200 lattices. In this case, we only 
require that the histogram of the random walk in the cor- 
responding energy segment is sufficiently flat without re- 
gard to the relative flatness over the entire energy range. 
In Fig. ||(b), the results for large lattices show clear dou- 
ble peaks for the canonical distributions at temperatures 
T^{L) = 0.70127 for L = 150 and T^{L) = 0.701243 for 
L = 200. The exact result is T^ = 0.701232.... for the 
infinite system. Considering the valley which we find for 
L = 200 is as deep as 9 x 10"^", we can understand why 
it is impossible for conventional Monte Carlo algorithms 
to overcome the tunneling barrier with available compu- 
tational resources. 

If we compare the histogram for L ~ 60 in Fig. ||(a) 
with that for L — 200 in Fig. ||(b), we see very clearly 
that the simulation effort for L — 200 (9.8 x 10^ visits per 
energy level) is even less then the effort for L = 60 (1.6 x 
10^ visits per energy level.). It is more efficient to perform 
random walks in relatively small energy segments than a 
single random walk over all energies. The reason is very 
simple, the random walk is a local walk, which means 
for a given Ei, the energy level for the next step only 
can be one of 9 levels in the energy range [Ei — A,Ei + 
4]. (for the Potts model discussed in this section). The 
algorithm itself only requires that the histogram on such 
local transitions is flat. (A single random walk, subject to 
the requirement of a flat histogram for all energy levels, 
will take quite long.) For random walks in small energy 
segments, we should be very careful to make sure that 
all spin configurations with energies in the desired range 
can be equally accessed so we restart the random walk 
periodically from independent spin configurations. 

An important question that must be addressed is the 
ultimate accuracy of the algorithm. One simple check is 
to estimate the transition temperature of the 2D Q = 10 
Potts model for L = oo since the exact solution is known. 
According to finite-size scaling theory, the "effective" 
transition temperature for finite systems behaves as: 



T,{L) = Te(cx3) + —^ 



(6) 



where Tc{L) and Tc(oo) are the transition temperatures 
for finite- and infinite-size systems, respectively, L is the 
linear size of the system and d is dimension of the lattice. 



In Fig. ^, the transition temperature is plotted as a 
function of L^"^. The data in the main portion of the 
figure are obtained from small systems(L = 10 ^ 30), 
and the error bars are estimated by results from multi- 
ple independent runs. Clearly the transition tempera- 
ture extrapolated from our simulational data is Tc(cxd) = 
0.7014 ± 0.0004 which is consistent with the exact solu- 
tion (Tc = 0.701232....) for the infinite system. To get 
an even more accurate estimate, and also test the accu- 
racy of the density of states from single runs for large 
systems, we performed single, long random walk on large 
lattices(L = 50 ~ 200). The results, plotted as a func- 
tion of lattice size in the inset of the figure, show that 
the transition temperature extrapolated from the finite 
systems is Tc(oo) = 0.701236 ± 0.000025 which is still 
consistent with the exact solution. 

We also compare our simulational result for the Q — 10 
Potts model with the existing numerical data such as 
estimates of transition temperatures and double peak 
locations obtained with the multicanonical simulational 
method by Berg and Neuhaus [gj and the Multibondic 
cluster algorithm by Janke and Kappler 0. All results 
are shown in Table 1. With our random walk simula- 
tional algorithm, we can calculate the density of states 
up to 200 X 200 within 10^ visits per energy level to ob- 
tain a good estimate of the transition temperature and 
locations of the double peaks. Using the multicanonical 
method and a finite scaling guess for the density of states, 
Berg et al. only obtained results for lattices as large as 
100 X 100 1^, and multibondic cluster algorithm data Q 
were not given for systems larger than 50 x 50. 

In section IV, the accuracy of our algorithm will be 
further tested by comparing thermodynamic quantities 
obtained for 2-dim Ising model with exact solutions. 



B. Thermodynamic properties of the Q 
model 



10 Potts 



One of advantages of our method is that the density 
of states does not depend on temperature; indeed with 
the density of states, we can calculate thermodynamic 
quantities at any temperature. For example, the internal 
energy can be calculated by: 



j:Eg{E)c 



-pE 



U{T) 



Y.9{E)e-PE 

E 



^{Eh 



(7) 



To study the behavior of the internal energy near Textc 
more carefully, we calculate the internal energy for L = 
60, 100 and 200 near Tc as presented in Fig. |^(a). A very 
sharp "jump" in the internal energy at transition tem- 
perature Tc is visible, and the magnitude of this jump is 
equal to the latent heat for the (first-order) phase tran- 
sition. Such behavior is related to the double peak dis- 
tribution of the first-order phase transition. When T is 



slightly away from Tc, one of the double peaks increases 
dramatically in magnitude and the other decreases. 

Since we only perform simulations on finite lattices, 
and use a continuum function to calculate thermody- 
namic quantities, all our quantities for finite-size systems 
will appear to be continuous if we use a very small scale. 
In the inset of Fig. ||(a) we use the same density of states 
again to calculate the internal energy for temperatures 
very close to Tc. On this scale the "discontinuity" at 
the first-order phase transition disappears and a smooth 
curve can be seen instead of a sharp "jump" in the main 
portion of Fig. ||(a). The discontinuity in Fig. ^(a) is 
simply due to the coarse scale, but when the system size 
goes to infinity, the discontinuity will be real. 

From the density of states we can also estimate the 
specific heat from the fluctuations in the internal energy: 



C(T) = 



dUjT) _ {E^)t - jE), 
dT 



rp2 



(8) 



In Fig. ^(b), the specific heat so obtained is shown as a 
fmiction of temperature. We calculate the specific heat in 
the vicinity of the transition temperature Tc. The finite- 
size dependence of the specific heat is clearly evident. We 
find that specific heat has a finite maximum value for a 
given lattice size L that, according to finite-size scaling 
theory for first-order transitions should vary as: 



c(T,T)T-'^(x/((T-Tc(oo))L'^) 



(9) 



where c{L, T) = C{L, T)/N is the specific heat per lattice 
site, L is the linear lattice size, d = 2 is the dimension of 
the lattice. T{L = cxj) = 0.70123.... is the exact solution 
for the Q = 10 Potts model ||]. In the inset of Fig. |(b), 
our simulational data for systems with L = 60, 100 and 
200 can be well fitted by a single scaling function, more- 
over this function is completely consistent with the one 
obtained from lattice sizes from L = 18 to T == 50 by 



standard Monte Carlo |36|. 

With the density of states, we not only can calcu- 
late most thermodynamic quantities for all temperatures 
without multiple simulations but can also access some 
quantities, such as the Gibbs free energy and entropy, 
which are not directly available from conventional Monte 
Carlo simulations. The free energy is calculated using 



Z= Y. e 

{ configurations} 

T = -A:Tlog(Z) 



-PE _ 



= J2g{E)e-^^ 



(10) 



Our results for the Gibbs free energy per lattice site is 
shown in Fig. ||(c) as a function of temperature. Since 
the transition is first-order the free energy appears to 
have a "discontinuity" in the first derivative at Tc. This 
is typical behavior for a first-order phase transition, and 
even with the fine scale used in the inset of Fig. ^(c), this 
property is still apparent even though the system is finite. 
The transition temperature Tc is determined by the point 



where the first derivative appears to be discontinuous. 
With a coarse temperature scale we can not distinguish 
the finite-size behavior of our model; however, we can 
see a very clear size dependence when we view the free 
energy on a very fine scale as in the inset of Fig. 0(c). 

The entropy is another very important thermodynamic 
quantity that cannot be calculated directly in conven- 
tional Monte Carlo simulations. It can be estimated by 
integrating over other thermodynamic quantities, such as 
specific heat, but the result is not always reliable since 
the specific heat itself is not easy to determine accu- 
rately, particularly considering the "divergence" at the 
first-order transition. With an accurate density of states 
estimated by our method, we already know the Gibbs 
free energy and internal energy for the system, so the 
entropy can be calculated easily: 



5(T) 



U{T) - F{T) 
T 



(11) 



It is very clear that the entropy is very small at low 
temperature and at T = is given by the density of states 
for the ground state. We show the entropy as a function 
of temperature in a wide region in Fig. ^.(d) 

The entropy has a very sharp "jump" at Tc, just as does 
the internal energy and such behavior can be seen very 
clearly in the inset of Fig. §(d), when we re-calculate the 
entropy near Tc . The change of the entropy at Tc shown 
in the figure can be obtained by the latent heat divided 
by the transition temperature, and the latent heat can be 
obtained by the jump in internal energy at Tc in Fig. ^(a). 

With the histogram method proposed by Ferrenberg 
and Swendsen [^, it is possible to use simulational data 
at specific temperatures to obtain complete thermody- 
namic information near, or between, those temperatures. 
Unfortunately it is usually quite hard to get accurate in- 
formation in the region far away from the simulated tem- 
perature due to difficulties in obtaining good statistics, 
especially for large systems where the canonical distribu- 
tions are very narrow. With the algorithm proposed in 
this paper, the histogram is "flat" for the random walk 
and we always have essentially the same statistics for all 
energy levels. Since the output of our simulation is the 
density of states, which does not depend on the temper- 
ature at all, we can then calculate most thermodynamic 
quantities at any temperature without repeating the sim- 
ulation. We also believe the algorithm is especially useful 
for obtaining thermodynamic information at low temper- 
ature or at the transition temperature for the systems 
where the conventional Monte Carlo algorithm is not so 
efficient. 



C. The tunneling time for the Q = 10 Potts model at 

To study the efficiency of our algorithm, we measure 
the tunneling time t, defined as the average number of 



sweeps needed to travel from one peak to the other and 
return to the starting peak in energy space. Since the his- 
togram that our random walk produces is flat in energy 
space, we expect the tunneling time will be the same 
as for the ideal case of a simple random walk in real 
space, i.e. t{Ne) ^ Ne, where Ne is the total number 
of energy levels. To compare our simulational results to 
those for the ideal case, we also perform a random walk 
in real space. We always use a fixed g{xi) = 1 in one- 
dimensional real space, where Xi is a discrete coordinate 
of position that can be chosen simply as 1, 2, 3 ,4 ...Ne- 
The random walk is a local random walk with transition 
probability p{xi -^ Xj) = 1/2, where Xj = Xi±i. We 
use the same definition of the tunneling time to measure 
the behavior of this quantity. The tunneling time for the 
ideal case satisfies the simple power law as t{Ne) ^ Ne" 
and the exponent a is equal to 1. (t is defined using the 
unit of sweep of Ne sites.) Our simulational data for ran- 
dom walks in energy space yield a tunneling time that is 
well described by the power law t ~ N" as shown in 
Fig. ^(a). The solid lines in the graph have the simple 
power law as t{L) ^ Ne , and we see that our simu- 
lation result is very close to the ideal case. Since our 
method needs an extra effort to update the density of 
states to produce a flat histogram during the random 
walk in energy space, the tunneling time is much longer 
than the real space case. Also because the tunneling time 
depends on the accuracy of the density of states, which 
is constantly modified during the random walk in energy 
space, it is not a well defined quantity in our algorithm. 
The tunneling time shown in Fig. Ha) is the overall tun- 
neling time which includes all iterations with the mod- 
ification factors from /o = C""^ ~ 2.71828... to the final 
modification factor /final = exp(10-**) ~ 1.00000001. 

We should point out that the two processes are not 
exactly the same, since the random walk in real space 
uses the exact density of sta.tes{g{xi) = 1). However 
the random walk in energy space requires knowledge of 
the density of states, which is a priori unknown. The 
algorithm we propose in this paper is a random walk with 
the density of states which is modified at each step during 
the walk in energy space. At the end of our random walk, 
the modification factor approaches 1, and the estimated 
density of states approaches the true value. The two 
processes are then almost identical. 

Conventional Monte Carlo algorithms (such as the 
heat-bath algorithm) have an exponentially fast grow- 
ing tunneling time. According to Berg's study in refer- 
ence H, the tunneling time obeys the exponential law 
t{L) = 1.46L2i5e°°s''^. The Multicanonical simula- 
tional method has reduced the tunneling time from an 
exponential law to a power law as t{Ne) ^ Ne"'. How- 
ever, the exponent a is as large as ~ 1.33 |q], which is 
far away from the ideal case a = 1. Only very recently, 
Janke and Kappler introduced the multibondic cluster 
algorithm, the exponent a is reduced to as small as 1.05 
for 2-dim ten state Potts model |Q. In Fig. ||(b), we 
show our result with the same quantities obtained with 



multicanonical method and heat bath algorithm in refer- 
ence g. We should point out that just like multicanon- 
ical simulational method, our algorithm has a power in- 
creasing tunneling time with a smaller exponent a. For 
small systems, our algorithm offers less advantage be- 
cause of the effort needed to modify the density of states 
during the random walk. Very recently, Neuhaus has 
generalized this algorithm to estimate the canonical dis- 
tribution for r < Tc, in magnetization space for the Ising 
model [^ . He found that for small systems the exponent 
for CPU time versus L for our algorithm and Multicanon- 
ical Ensemble simulations are almost identical. However 
we see that our algorithm is very efficient for large sys- 
tems, especially for L > 100. Our results in Fig. Mb) are 
only for single range random walks, and multiple range 
random walks have been proven more efficient for larger 
systems. 



IV. APPLICATION TO A SECOND ORDER 
PHASE TRANSITION 

The algorithm we proposed in this paper is very effi- 
cient for the study of any order phase transitions. Since 
our method is independent of temperature, it reduces 
the critical slowing down at the second-order phase tran- 
sition Tc and slow dynamics at low temperature. We 
estimate the density of states very accurately with a fiat 
histogram, the algorithm will be a very efficient for gen- 
eral simulational problems by avoiding the need for mul- 
tiple simulations at multiple temperatures. 

To check the accuracy and convergence of our method, 
we apply it to the 2D Ising model with nearest neighbor 
interactions on a L x L square lattice. This model pro- 
vides an ideal benchmark for new algorithms p,p8[ and 
is also an ideal laboratory for testing theory p0| , |39[ . This 
model can be solved exactly, therefore we can compare 
our simulational results with exact solutions. 

With the exact solution for the partition function on 
finite-size systems ||lO| , and expansion of the expression 
by Mathematica, the density of states for the 2D Ising 
model on a square lattice can be obtained exactly [[l()[ . 
Beale llOJ obtained the exact density of states up to i = 
32; and because of the memory and speed limitation of 
present computers, we could only extend the calculation 
up to L = 50 with the Mathematica program provided 
by Beale. With the algorithm proposed in this paper, the 
density of states for the Ising model is estimated for L — 
50. The final modification factor for our random walk 
was 1.000000001 for L = 50. In Fig. |(a) the simulational 
densities of states are compared with exact results. With 
the logarithmic scale used in the figure simulational data 
and exact solution overlap perfectly with each other for 
L = 50. In the inset of the figure, we show the relative 
error e, which is generally defined by the ratio between 
error of the simulational data and exact values for any 
quantity X: 



eiX) 



\X, 



x« 



X, 



(12) 



e{\og{g)) for most of the region is smaller than 0.1%. 
Such errors for low energy levels are directly related to 
the errors for the thermodynamic quantities calculated 
from the density of states. The average relative error is 
0.019% for L = 50. It is possible to estimate the den- 
sity of states for small systems with the broad histogram 
method |21-23|. Recent broad histogram simulational 
data M] for the 2D Ising model on a 32 x 32 lattice 
with 10^ MC sweeps yielded an average deviation of the 
microcanonical entropy of about 0.08 % from the exact 
solution [M. With our algorithm we obtain an average 
error as small as 0.035 % on the 32 x 32 lattice with 7x 10^ 
MC sweeps. Procedures for allowing / ^ 1 have been ex- 
amined by Hiiller pd[ | who used data from two densities 
of states for two different values of / to extrapolate to 
/ = 1. However, his data for a small Ising system yield 
larger errors than our direct approach. The applicability 
of his method to large systems also needs a more detailed 
study. 

The absolute density of states in Fig. @(a) is obtained 
by the condition that the number of ground states is 2 
for the 2D Ising model (all up or down). This condition 
guarantees the accuracy of the density of states at low 
energy levels which are very important in the calculation 
of thermodynamic quantities at low temperature. With 
this condition, when T = 0, we can get exact solutions 
for internal energy, entropy and free energy when we cal- 
culate such quantities from the density of states. If we 
apply the condition that the total number of states is 2^ 
for the ferromagnetic Ising model, we can not guaran- 
tee the accuracy of the energy levels at or near ground 
states because the rescaled factor is dominated by the 
maximum density of states. 

In Fig. ^(b) , we show our estimation of the density of 
states of Ising model on 256 x 256 lattice. Since the 
density of states for E > has almost no contribu- 
tion to the canonical average at finite positive temper- 
ature, we only estimate the density of states in the re- 
gion E/N G [-2, 0.2] out of the whole energy [-2, 2]. To 
speed up our calculation, we divide the desired energy 
region [-2, 0.2] into 15 energy segments, and estimate 
the density of states for each segment with independent 
random walks. The modification factor changes from 
/o = el ~ 2.71828... to /fi„ai = exp(10-^) ~ 1.0000001.... 
The resultant density of states can be joined from adja- 
cent energy segments. To reduce the boundary effects of 
the random walk on each segments, we keep about sev- 
eral hundred overlapping energy levels for random walks 
on two adjacent energy segments. The histograms of ran- 
dom walks are shown in the inset of this figure. We only 
require a flat histogram for each energy segment. To re- 
duce the error of the density of states relevant to the 
accuracy of the thermodynamic quantities near T^ we 
optimize the parameter and perform additional multiple 
random walks for the energy range E/N E [—1.8,-1] 



with same number of processors. For this we use the 
density of states obtained from the first simulations as 
starting points and continue the random walk with mod- 
ification factors changing from exp(10~^) ~ 1.000001 to 
exp(10"^) ~ 1.000000001. The total computational ef- 
fort is about 9.2 x 10^ visits on each energy levels. Note 
that the total number of possible energy levels is iV — 1 
and we perform random walks only on [-2, 0.2] out of [-2, 
2]. The real simulational effort is about 6.1 x 10^ MC 
sweeps for the Ising model with L = 256. With the pro- 
gram we implemented, it took about 240 CPU hours on 
a single IBM SP Power3 processor. 

For L = 256, we perform multiple random walks on 
different energy ranges, and one problem arises, that 
is the error of the density of states due to the random 
walk in a restricted energy range. Since the exact den- 
sity of states for large systems is not available, we use 
L = 32 to study such effects. We perform three indepen- 
dent random walks in the ranges E/N = [—1.7,-1.2], 
E/N = [-1.8,-1.1] and E/N = [-1.9,-1.0] to calcu- 
late the densities of states on these ranges. In Fig. 0(a), 
we show the errors of our simulation results from the ex- 
act values. We make our simulational densities of states 
match up with the exact results at the left edges. It 
is very clear that the width of the energy range of the 
random walks is almost not relevant to the errors of the 
density of states. The reason is that the random walks 
only require the local histogram to be flat as we discussed 
in the previous section. But the errors for the last two 
densities on the right edges are significantly larger than 
the others. The problem is the boundary effect because 
in this simulation we treat the densities of states at edges 
as same as those away from edges. That is not correct. 
Another reason is due to the flat histogram and high den- 
sity at right side edges. Our simulational effort for the 
densities at the right edges is not enough compared to 
those at the left side. This is also the reason why we 
have not seen big errors at left edges. 

To study the influence of the errors of the densities of 
states on the thermodynamic quantities calculated from 
them, in the energy range that we perform random walks, 
we replace the exact density of states with the simula- 
tional density of states. In Fig. 0(b), the specific heat 
calculated from such density of states is showed as a func- 
tion of temperature. We also show the exact value with 
the simulational data, the difference is obvious. 

To reduce the boundary effect, we delete the last two 
density entries, and insert them into the exact density 
of states again, then the results in Fig. M(c) are much 
better. 

With our test in the three different ranges of energy, it 
is quite safe to conclude that boundary effect will not be 
present in our multiply random walks if we have a couple 
of energy levels overlap for adjacent energy ranges. In 
our real simulations for large systems, we have hundreds 
of overlaping energy levels. 

Since the exact density of states is only available on 
small systems (L < 50), it is not so interesting to com- 



pare the simulational density of states itself. The most 
important thing is the accuracy of estimations of ther- 
modynamic quantities calculated from such density of 
states on large systems. With the density of states on 
large systems, we apply canonical average formulas to 
calculate internal energy, specific heat, Gibbs free energy 
and entropy. Ferdinand and Fisher l^^l obtained the ex- 
act solutions of above quantities for 2D Ising model on 
finite-size lattices. Our simulational results on finite-size 
lattice can be compared with those exact solutions. 

In Fig. @(a), we show the internal energy as a function 
of temperature. The internal energy is estimated from 
the canonical average over energy of the system as the 
equation (^. We also draw the exact solution in the same 
graph. The exact and simulational data perfectly overlap 
with each other in a wide temperature region from T = 
to T — 8. Since no difference is visible with the scale 
used in Fig. 0(a) , a more stringent test of the accuracy 
is provided by the inset which shows the relative errors 
s{U). With the density of states obtained with our algo- 
rithm, the relative error of simulational internal energy 
for L = 256 is smaller than 0.09% for the temperature re- 
gion from r = to 8. From eqn. (^ it is very clear that 
the canonical distribution serves as a weighting factor, 
and since the distribution is very narrow, U(T) is only 
determined by a small portion of the density of states. 
(For the i = 50 2D Ising model at Tc, only the density 
of states for E/N e [—1.6,-1.2] contributes in a major 
way to the calculation.) Therefore the error e{U) is also 
determined by the errors of the density of states in the 
same narrow energy range. 

With the density of states and equation (|8[), we also can 
calculate the specific heat per lattice site as a function of 
temperature. Both simulational data and exact results 
near T^ are shown in Fig. |(b) for L = 64, 128 and 256 
Ising model in the vicinity of Tc. Within the resolution of 
the figure, we can see the difference between the simula- 
tional data and exact solutions. In the inset of Fig. H(b), 
we present relative errors for our simulational data as a 
function of temperature for L = 256. The errors for the 
specific heat for the Ising model on a 256 x 256 lattice are 
smaller than 4.5% in all temperature region T < 8. Very 
recently, Wang, Tay and Swendsen |38) estimated the 
specific heat of the same model on a 64 x 64 lattice by the 
transition matrix Monte Carlo re- weighting method [ p2[ , 
and for a simulation with 2.5 x 10^ MC sweeps, the max- 
imum error in temperature region T e [0, 8] was about 
1%. When we apply our algorithm to the same model 
on the 64 x 64 lattice, with a final modification factor of 
1.000000001 and a total of 2 x 10^ MC sweeps on sin- 
gle processor, the errors in the specific heat are reduced 
below 0.7% for all temperatures. The relatively large er- 
rors at low temperature reflect the small values for the 
specific heat at low temperature. According to the sim- 
ulational data by broad histogram method, the errors in 
specific heat are very large even for systems as small as 
32 X 32 pl|-p3[. Only very recently, they have reduced 
the error near Tc to a small value for L = 32 [Es. 



In Fig. ||(c), the free energy for L = 256 is plotted 
as a function of temperature. For comparison the exact 
solution En] is shown in the same figure. As expected, 
simulational and exact data overlap perfectly within the 
resolution of the figure. Since the system has a second- 
order phase transition, unlike the Q = 10 Potts model, 
the first derivative of the free energy is a continuous func- 
tion of temperature. The result very close to Tc{L — 256) 
is shown in the inset of Fig. H(c) for L = 256. The behav- 
ior of the free energy near Tc is quite different from that 
of the first-order phase transition in Fig. 0(c). According 
to our calculation, the relative errors for all tempera- 
tures from T = to 8 are smaller than 0.0008%, which 
means that our simulational data agree almost perfectly 
with the exact solution [g^. Since we use the condition 
that the number of ground states is g{E = —27V) = 2 to 
normalize the density of states, the errors at low temper- 
atures are extremely small. 

The entropy of the 2D Ising model can be calculated 
with the equation(pT|). In Fig. ^(d), the simulational data 
and exact results are presented in the same figure. With 
the scale in the figure, the difference between our simu- 
lational data and exact solutions are not visible. In the 
inset of Fig. @(d), the relative errors of our simulational 
data are plotted as a function of temperature. For the 
Ising model on a 256 x 256 lattice, the relative errors 
are smaller than 1.2% for all temperature range. We 
also notice that the errors near T — decrease dramati- 
cally. The reason is the condition we use to normalize the 
density of states. Very recently, with the flat histogram 
method |2|] and the broad histogram method [pll-p3[, 
the entropy was estimated with 10^ MC sweeps for the 
same model on 32 x 32 lattice; however, the errors in 
reference |24] are even much bigger than our errors for 
256 X 256! 



V. APPLICATION TO 3D ±J EA MODEL 

Spin glasses |44| are magnetic systems in which the 
interactions between the magnetic moments produce 
frustration because of some structural disorder. One 
of the simplest theoretical models for such systems is 
the Edwards- Anderson model U^ (EA model) proposed 
twenty five years ago. For such disordered systems, an- 
alytical methods can provide only very limited informa- 
tion, so computer simulations play a particularly impor- 
tant role. However, because of the rough energy land- 
scape of such disordered systems, the relaxation times of 
the conventional Monte Carlo simulations are very long. 
The dynamical critical exponent was estimated as large 
as z ~ 6 |46| - |48[ . Normally simulations can be performed 
only on rather small systems, and many properties con- 
cerning the spin glasses are still left unclarified |4S|-[5q]. 

In this paper, we consider the three-dimensional ±J 
Ising spin glass EA model. The model is defined by the 
Hamiltonian 
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where a is an Ising spin and the couphng J^ is quenched 
to ±1 randomly. The summation runs over the nearest- 
neighbors {i,j) on a simple cubic lattice. 

One of most important issues for a spin glass model is 
the low temperature behavior. Because of the slow dy- 
namics and rough phase space landscape of this model, 
it is also one of most difficult problems in simulational 
physics. The algorithm proposed here is not only very 
efficient in estimating the density of states but also very 
aggressive in finding the ground states. From a random 
walk in energy space, we can estimate the ground state 
energy and the density of states very easily. For a spin 
glass system, after we finish the random walk, we can ob- 
tain the absolute density of states by the condition that 
total number of states is 2^. The entropy at zero tem- 
perature can be calculated from either ^o = ln{g{Eo)) or 
lim ^~'^ , where Eq is the energy at ground states. Both 

relations will give the same result since U and F are cal- 
culated from the same density of states. Our estimates 
for So = Sq/N and cq — Eq/N per lattice site, listed 
in Table 2, agree with the corresponding estimates made 
with the multicanonical method. With our algorithm, 
we can estimate the density of states up to L = 20 by a 
random walk in energy space for few hours on a 400MHz 
processor. 

If we are only interested in the quantities directly re- 
lated to the energy, such as free energy, entropy, internal 
energy and specific heat, one dimensional random walk 
in energy space will allow us to calculate these quan- 
tities with a high accuracy as we did in the 2D Ising 
model. However for spin glass systems, one of the most 
important quantities is the order parameter which can 
be defined by pi 
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(T) = hm hm q{T,t), q{T,t) = (V (T,(0)a,(t)/iV). 



i^l 



(14) 



Here, iV = L^^ is the total number of the spins in the 
system, L is the linear size of the system, q{T,t) is the 
auto-correlation function, which depends on the temper- 
ature T and the evolution time t, and q{T, 0) = 1. When 
t — + oo, q{T,t) becomes the order parameter of the spin 
glass. This parameter takes the following values 



-EA 



= 1 if T = 
{T){ = if T > Tg 
T^Oif 0<T<rj 



(15) 



The value at T = can be different from 1 in the case 
where the ground state is highly degenerate. 

In our simulation, there is no temperature introduced 
during the random walk. And it is more efficient to per- 
form a random walk in single system than two replicas. 
So the order-parameter can be defined 



q^{J2a^a,/N). 



(16) 



i=l 



where {cr^} is one of spin configurations at ground states 
and {<Ji} is any configuration during the random walk. 
The behavior of q we defined above is basically the same 
as the order-parameter defined by the Edwards and An- 
derson p5| . It is not exact same order-parameter defined 
by Edwards and Anderson, but was used in the early 
numerical simulations by Binder et al. [ p7| , p8[ . 

After first generating a bond configuration we perform 
a one-dimensional random walk in energy space to find a 
spin configuration {a'^} for the ground states. Since the 
order-parameter is not directly related to the energy, to 
get a good estimate of this quantity we have to perform 
a two-dimensional random walk to obtain the density of 
states G{E, q) with a flat histogram in E-q space. This 
also allows us to overcome the barriers in parameter space 
(or configuration space) for such a complex system. The 
rule for the 2D random walk is the same as ID random 
walk in the energy space. 

In Fig. g, we show the histogram of the 2D random 
walk in energy-order parameter space, which is very flat. 
With the density of states G{E, q), we can calculate any 
quantities as we did in the previous sections. It is very 
interesting to study the roughness of this model. First 
we study the canonical distribution as a function of the 
order-parameter : 



P(g,r) = ^G(i?,g)f 



-E/kuT 



(17) 



In Fig. &Q(a), we show a 3D plot for the canonical 
distribution at different temperatures for one bond con- 
figuration of L = 6 EA model. At low temperatures, 
there are four peaks, and the depth of the valleys be- 
tween peaks depends upon temperature. When the tem- 
perature is high, the multiple peaks converge to a single 
central peak. Because we use the linear scale to show our 
result in the Fig. EO(a) . It is not clear how deep the dips 
among peaks are. In Fig. |lO|(b), we show the canonical 
distribution using logarithmic scale for the same distri- 
bution but only at T = 0.5, and we find that the dips 
are as deep as 10"'*. We also noted there actually are six 
peaks, but the plot with linear scale does not show all 
of them because two are as small as 10"'^ compared to 
other four peaks. 

In Fig. |ll|(a), we show the roughness of the canoni- 
cal distribution for another realization on an 8'^ lattice. 
Because of the wide variation in the distribution at low 
temperature we used a logarithmic scale: the relative size 
of dips are as deep as lO^'^" at T = 0.1. There are several 
local minima even at high temperatures. With conven- 
tional Monte Carlo simulations, it is almost impossible 
to overcome the barriers at the low temperature, so the 
simulation will get trapped in one of the local minima as 
shown in the figure. With our algorithm all states will 
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be visited with more or less the same probabiHty and 
trapping is not a problem. 

With the density of states G{E, q), we also can calcu- 
late the energy landscape by: 



U{q,T) 



j:EG{E,q)e-^^ 
j:GiE,q)e-PE 

E.q 



(18) 



In Fig. 0(b), we show the internal energy as a function 
of order-parameter for temperatures T = 0.1 ^ 2.0. We 
find that the landscape is a very rough at low temper- 
atures. The roughness of the energy landscape agrees 
with the one for canonical distribution. But the maxima 
in energy landscape are corresponding to the minima ap- 
proximately in the canonical distribution. 

As we already noted in the previous paragraph, the 
roughness of the landscape of the spin glass model makes 
the conventional Monte Carlo simulation extremely dif- 
ficult to apply. Therefore, even a quarter of a cen- 
tury after the model was proposed, we even can not 
conclude whether there is a finite phase transition be- 
tween glass phase and disordered phase. With Monte 
Carlo simulations on a large system (64^ x 128) and a 
finite-size scaling analysis on a small lattice, Marinari et 
al. |59| expressed doubt about the existence of the "well- 
established" finite-temperature phase transition of the 
3D Ising spin glass |^,Q. Their simulational data can 
be described equally well by a finite-temperature tran- 
sition or by a T = singularity of an unusual type. 
Kawashima and Young's simulational data could also 
not rule out the possibility of Tg = [g0|. Thus even 
the existence of the finite-temperature phase transition 
is still controversial, and thus the nature of the spin glass 
state is uncertain. Although the best available computer 
simulation results [15 54 60] have been interpreted as a 



mean-field like behavior with replica-symmetry break- 
ing(RSB) 1^, Moore et al. showed evidence for the 
droplet picture |£2| of spin glasses within the Migdal- 
Kadanoff approximation. They argued that the failure 
to see droplet model behavior in Monte Carlo simula- 
tions was due to the fact that all existing simulations are 
done at temperatures too close to transition temperature 
so that system sizes larger than the correlation length 
were not used. As discussed in the previous paragraph, 
the lower the temperature is, the rougher the canoni- 
cal distribution and energy landscape are; hence, it is 
almost impossible for conventional Monte Carlo meth- 
ods to overcome the barrier between local minima and 
globe minima. It is possible to heat the system up to in- 
crease the possibility of escape from local minima by sim- 
ulated annealing and the more recent simulated temper- 
ing method p3[ and parallel tempering method p3,R5[, 
but it is still very difficult to perform equilibrium simu- 
lations at low temperatures. Very recently Hatano and 
Gubernatis proposed a bivariate multicanonical Monte 
Carlo method for the 3D ± J spin glass model, and their 
result also favors the droplet picture liaJHfl]. Marinari, 



Parisi et al. argued, however, that the data were not 
thermalized |60|. The nature of spin glasses thus remains 
controversial |m] . 

The algorithm proposed in this paper provides an al- 
ternative for the study of complex systems. Because we 
need to calculate the order-parameter with high accuracy, 
and this quantity is not directly related to the energy, we 
need to perform a random walk in the two-dimensional 
energy - order parameter space. After we estimate the 
density of states in this 2D space, we can calculate the 
order-parameter at any temperature from the canonical 
average. In |2|(a), we show our results for the iD EA 
model for L = 4, 6 and 8. Because we need to perform 
a 2D random walk with a total of about L^ states, the 
simulation is only a practical for a small system(L < 8). 
The results in the figure are the average over 100 real- 
izations for L = 4, 50 realizations for L = 6 and 20 for 
L = 8. 

We notice that the behavior of {q{T)) is very similar 
to the magnetization (the order-parameter for the Ising 
model), but the finite value at low temperature is not 
necessarily equal to 1 because of the high degeneracy of 
the ground state for the spin glass model. The fluctuation 
of the order-parameter at the different temperatures for 
L = 4, 6 and 8 is shown in the inset of the figure. 

To estimate the transition temperature of the spin 
glass system, we calculated the fourth order cumulant 
as a function of temperature. In Fig. |l^(b), we show our 
simulational results for L = 4, 6 and 8. All curves clearly 
cross around Tg = 1.2. Below this temperature, the 
spin configurations are frozen into some disorder ground 
states and the order parameter assumes a finite value. 
Above this temperature Tg, the system is in a disordered 
states and the order parameter vanishes. 

One complication for simulation of such random sys- 
tems is the determination of the relative importance of 
the error due to the simulation algorithm and the error 
due to the finite sampling of bond distributions. From 
Fig. |l2| (a) and (b), we can not tell what the origin of the 
error bars is so we also performed multiple independent 
simulations for the same bond configuration on a L = 6 
3D EA model. We found that the statistical errors for 
the order parameter and the fourth order cumulant from 
these simulations were much smaller than the error bars 
shown in Fig. |2| (a) and (b) for all temperatures. We 
conclude that the error bars in the figure arise almost 
completely from the randomness of the system. 

The computational resources devoted here to the EA 
model were not immense. All our simulations for one 
bond configuration (L = 4, 6 and 8) were performed 
within two days on (multiple) Linux machines (200 — 
800MHz) in the Center for Simulational Physics. This 
effort should thus be viewed as a feasibility study, and 
substantially more effort would be required to determine 
the nature of the spin glass phase or to estimate the tran- 
sition temperature with high accuracy. Nonetheless, we 
believe that these results show the applicability of our 
method to systems with a rough landscape. Because of 
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the number of states is about 7V^ for 2D random walks, 
such calculations not only require huge memory during 
the simulation but also substantial disk space to store 
the density of states for the later calculation of thermo- 
dynamic quantities. 



VI. DISCUSSION AND CONCLUSION 

In this paper, we proposed an efficient algorithm to cal- 
culate the density of states directly for large systems. By 
modifying the estimate at each step of the random walk 
in energy space and carefully controlling the modification 
factor, we can determine the density of states very accu- 
rately. Using the density of states, we can then calculate 
thermodynamic quantities at essentially any temperature 
by applying simple statistical physics formulas. An im- 
portant advantage of this approach is that we can also 
calculate the Gibbs free energy and entropy, quantities 
that are not directly available from conventional Monte 
Carlo simulations. 

We applied our method to the 2D Q = 10 Potts model 
which demonstrates a typical first-order phase transition. 
By estimating the density of states with lattices as large 
as 200 X 200, we calculated the internal energy, specific 
heat, free energy and entropy in a wide temperature re- 
gion. We found a typical first-order phase transition with 
a "discontinuity" for the internal energy and entropy at 
Tc. The first derivative of the free energy also shows such 
a discontinuity at T^. The transition temperature esti- 
mated from simulational data is consistent with the exact 
solution. 

We also applied our algorithm to the 2D Ising model, 
which shows a second-order phase transition. The den- 
sity of states obtained by the end of our simulations was 
compared directly with the exact solution on 50 x 50 lat- 
tice. The relative errors for most important energy levels 
are less than 0.019%. It was also possible to calculate 
the density of states for a 256 x 256 lattice with a com- 
putational effort of 6.1 x 10^ Monte Carlo sweeps. With 
the accurate density of states, we calculated the inter- 
nal energy, specific heat, Gibbs free energy and entropy. 
For all temperatures between T = and T = 8, the rel- 
ative errors are smaller than 0.09% for internal energy, 
0.0008% for free energy, 1.2% for entropy and 4.5% for 
specific heat. 

We should point out that our simulational results for 
L = 256 are close to the same accuracy as the specific 
heat of the 2D Ising model for L — 64 with the transition 
matrix Monte Carlo re- weighting method with 2.5 x 10^ 
MC sweeps pq] . Our estimate of the entropy for the 2D 
Ising model for the L — 256 lattice is even more accurate 
than the results obtained for the same quantity for L — 
32 with the broad histogram method and flat histogram 
method [e3, which needed about 10^ MC sweeps. Our 
simulational effort for L = 256 was 6.1 x 10^ MC sweeps. 

The algorithm was also applied with success to the 3D 



±J EA spin glass model for which we could determine 
the roughness of the energy landscape and canonical dis- 
tribution in order-parameter space. The internal energy 
and entropy at zero temperature were estimated up to a 
lattice size 20"^, and the transition temperature was esti- 
mated as about Tg = 1.2. 

In this paper, we only concentrated the random walk 
in energy space (and order-parameter space); however, 
the idea is very general and we can apply this algorithm 
to any parameters [|l^ . The energy levels of the models 
treated here are perfectly discrete and the total num- 
ber of possible energy levels is known before simulation, 
but in a general model such information is not available. 
Since the histogram of the random walk with our algo- 
rithm tends to be flat, it is very easy to probe all possible 
energies and monitor the histogram entry at each energy 
level. For some models where all possible energy levels 
can not be fitted in the computer memory or the energy 
is continuous, e.g. the Heisenberg model, we may need 
to discretize the energy levels. According our experience 
on discrete and continuous models, if the total number 
of possible energies is around the number of lattice sites 
N, the algorithm is very efficient for studying both first- 
or second-order phase transitions. 

In this paper, we only applied our algorithm to sim- 
ple models, but since the algorithm is very efficient even 
for large systems it should be very useful in the studies 
of general, complex systems with rough landscapes. It 
is clear, however, that more investigation is needed to 
better determine under which circumstances our method 
offers substantial advantage over other approaches and 
we wish to encourage the application of this approach to 
other models. 
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FIG. 1. Density of states g{E) for the 2-dim Q = 10 Potts 
model as a function of energy per lattice site E/N . With 
the scale in the figure, the errors of the simulational data are 
within the width of the lines. 

Figl 



40000 



(b) 



Fig 2(b) 



30000 
? 20000 



10000 




10" 



10' 



K 10- 

CD 

lD 
'(D 
g 10 ' 

10" 



10" 





l'\ 


A 


10* 

X 

10* 
10' 

1 ^10= 

; X 

10* 






'\ 


(b) \ 

\ '1' 


U150 


\\ui5o// 
\U200 / 


2 -1.5 -1 -0.5 


V 

L=200 






2 -1.5 -1 -0.5 

E/N 



-1.5 



-1 
E/N 



-0.5 



FIG. 2. The canonical distributions at the transition tem- 
perature P{E,Tc) = p(£;)e-^/^^'^'= for the Q = 10 Potts 
model for (a) L = 60, 80 and 100 (single random walk) and 
(b) L = 150 and 200 (multiple random walks). The insets 
show the histograms of the random walks to estimate the 
densities of states. 
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FIG. 3. Extrapolation of finite lattice "transition temper- 
atures" for the Q = 10 Potts model. 
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FIG. 4. Thermodynamic quantities calculated from the 
density of states for the Q = 10 Potts model: (a) internal 
energy, (b) specific heat and the finite-size scaling function, 
(c) Gibbs free energy and (d) entropy 
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FIG. 5. Tunneling times r for the Q = 10 Potts model: 
(a) comparison of r for our random walk algorithm in en- 
ergy space and for an ideal random walk in real space; (b) 
comparison of the tunneling times for our algorithm, for the 
multicanonical ensemble method and for the heat bath algo- 
rithm. The solid lines show the ideal case with t{Ne) ~ Ne- 
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FIG. 6. Density of states {log{g{E))) of the 2D Ising model 
for (a) 1/ = 50 (single random walks) and (b) L — 256 (mul- 
tiple random walks). The relative errors of the simulational 
densities of states are shown in the inset of (a) . The overall 
histogram of the random walk is shown in the inset of (b). 
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FIG. 7. Boundary effects for random walks in different en- 
ergy ranges for the 2D Ising model, (a) errors in the density of 
states; (b) specific heats calculated from the density of states; 
and (c) the specific heat if the two highest energy entries in 
the density of states are deleted. 
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FIG. 8. Thermodynamic quantities for the 2D Ising model 
calculated from the density of states. Relative errors with 
respect to the exact solutions by Ferdinand and Fisher are 
shown in the insets, (a) Internal energy, (b) specific heat, (c) 
Gibbs free energy and (d) entropy 
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FIG. 9. The histogram of the two-dimensional random 
walk in energy and order-parameter space for the 3D EA 
spin-glass model. 



H(E,q) 




FIG. 10. a) Overview of the rough topology of the canon- 
ical distribution in the order-parameter space for one bond 
configuration of the 3-D EA model on an L = 6 simple cubic 
lattice, b) The logarithmic plot for the canonical distribution 
as a function of the order-parameter for the 3D EA model at 
the temperature T = 0.5. 

(a) 



0.04 



Fig10(a) 



0.03 



S 0.02 




0.01 



(b) 



10" 



10' 



10' 



10" 



10"' 



10" 




20 



FIG. 11. (a) The canonical distribution P{q) as a function 
of the order-parameter for one bond configuration of the 3D 
EA model on an L=8 simple cubic lattice at the temperature 
T = 0.1 ~ 2.0. (b) Corresponding energy landscapes U{q,T) 
at different temperatures. 
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Fig 11 (a) 




FIG. 12. Properties of the 3-D EA spin glass model calcu- 
lated from the density of states G{E, q) resulting from a 2D 
random walk in energy - order parameter space: (a) The order 
parameter vs temperature; (b) The temperature dependence 
of the fourth order cumulant of the order parameter. The 
cumulants for different lattice sizes cross around Tc = 1.2. 

(a) 

Fig 12(a) 



0.8 



0.6 



0.4 



0.2 



^~X 



i; 





^^ 


15 




/' i 


10 




,.i-. '\ 


5 






^ 




^jy^ ^~~~^=i=a,=..^ 





^-^ 1 1 


) 1 2 3 

T 


■■=3=,. ^^^^^^^^^^-^ — ___ 


■"--, 








~^ 



2 

T 



(b) 



Fig 11(b) 



(b) 



Fig 12(b) 




21 



TABLE I. Estimates of "transition temperature" Tc and positions of double peaks _B™"'', _E™"^ for the Q = 10 Potts model 
with our method, Multicanonical (MUCA) ensemble[5] and Muhibondic (MUBO) cluster algorithm [6]. ET"'' and E^"^ are 
the energy per lattice site at the two peaks of canonical distribution at Tc . 



size 


Our method 


MUCA 


MUBO 


L 


Te 


rpmax 

El 


rpmax 
-t'2 


Tc 


rpmax 


rpmax 

E2 


Tc 


rpmax 


rpmax 


12 


0.70991 


0.8402 


1.7013 


0.710540 


0.806 


1.688 


0.7105402 


0.833 


1.72 


16 


0.70653 


0.8694 


1.6967 


0.706544 


0.844 


1.676 


0.7065144 


0.867 


1.71 


20 


0.70511 


0.8925 


1.6875 








0.7047891 


0.883 


1.69 


24 


0.70362 


0.8940 


1.6765 


0.703730 


0.908 


1.698 








26 


0.70317 


0.9002 


1.6805 








0.7034120 


0.908 


1.682 


30 


0.70289 


0.9233 


1.6888 














34 


0.70258 


0.9343 


1.6732 


0.702553 


0.927 


1.683 


0.7025530 


0.921 


1.676 


40 


0.70239 


0.9337 


1.6731 














50 


0.70177 


0.9416 


1.6776 








0.7018765 


0.940 


1.674 


60 


0.70171 


0.9522 


1.6733 














70 


0.70153 


0.9519 


1.6717 


0.701562 


0.9511 


1.670 








80 


0.70143 


0.9576 


1.6701 














90 


0.70141 


0.9551 


1.6727 














100 


0.70135 


0.9615 


1.6699 


0.701378 


0.9594 


1.6699 








120 


0.70131 


0.9803 


1.6543 














150 


0.70127 


0.9674 


1.6738 














200 


0.70124 


0.9647 


1.6710 














oo 


0.701236 ± 0.000025 






exact 


0.701232... 







TABLE II. Estimates of entropy (sq) and internal energy 
our method, multicanonical Method(MUCA) [17]. 



eo) per lattice site at zero temperature for the 3D EA model by 



size 


Our Method 




MUCA 




L 


So 


eo 


So 


eo 


4 


0.075 ± 0.027 


-1.734 ±0.006 


0.0724 ± 0.0047 


-1.7403 ±0.0114 


6 


0.061 ± 0.025 


-1.767 ±0.024 


0.0489 ± 0.0049 


-1.7741 ±0.0074 


8 


0.0493 ± 0.0069 


-1.779 ±0.016 


0.0459 ± 0.0030 


-1.7822 ±0.0081 


12 


0.0534 ± 0.0012 


-1.780 ±0.012 


0.0491 ± 0.0023 


-1.7843 ±0.0030 


16 


0.0575 ± 0.0037 


-1.7758 ±0.0041 






20 


0.0556 ± 0.0034 


-1.7745 ±0.0043 
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